Mapping Basics
Advanced Data Analytics

1 Mapping in R

To be fair, cartography and spatial analysis is not one of my favorite things to do in R. While nearly every statistical analysis you want to do has several packages that can accomplish the same task, rarely do they rely one drastically different types of data. This is why a lot of spatial analysis is still performed in programs such as ArcGIS or QGIS. However, moving data from one program to another often requires reformatting the data which can also be cumbersome. So maintaining the data within a single program to manage, format, analyze, and visualize helps to create a usable, repeatable workflow. That’s why data analysis and visualization in R is becoming standard practice in many disciplines.

1.1 Data Types

When using R for spatial analysis there are a number of graphics packages that can be used for everything from basic mapping to interactive maps. A few of these packages you may already be familiar with such as plot from base R, ggplot2, plotly, leaflet. Many of these packages are dependent on other packages that help you manage, organize, and format spatial data. Packages such as sf “Simple Features for R”, sp “Classes and Methods for Spatial Data”, and rgdal “Bindings for the ‘Geospatial’ Data Abstraction Library” provide access to spatial data and spatial data formats for a number of other packages. In order to reduce the overall complexity of this exercise, we will not discuss all of the data formats or packages that can be used to analyze spatial data here. However, I will provide links to useful online resources in the exercise README document.

1.2 Packages

For this exercise we will use a few common packages for working with spatial data. Once again, many of these packages will require dependencies so installing them may take some time.

ggsn | leaflet | mapdata | maptools | OpenStreetMap | rgdal | smoothr | sf | sp | tidyverse | tmap

In order to avoid forcing these packages to your computer, search in the Packages Tab of RStudio to see if you need to install them. Once you have installed them you can continue with the script below:

packages<-c("ggsn","leaflet","mapdata","maptools","OpenStreetMap","rgdal","smoothr","sf","sp","tidyverse","tmap")
sapply(packages, require, character.only=T)

         ggsn       leaflet       mapdata      maptools OpenStreetMap         rgdal       smoothr            sf 
         TRUE          TRUE          TRUE          TRUE          TRUE          TRUE          TRUE          TRUE 
           sp     tidyverse          tmap 
         TRUE          TRUE          TRUE 

These packages will provide the minimum functionality required to complete this exercise. As you experiment with various styles of mapping you may require additional packages to be installed at a later time.

2 Simple Maps

There are a number of very simple maps you can create that uses data included in the package downloads. For this first example we are going to create a map that depicts the location of Austin Peay State University. To begin we need to load some data from ggplot2 using the map_data command.

state <- map_data("state")
county <- map_data("county")
apsu_point <- data.frame("x" = -87.353069, "y" = 36.533654)

This quick script loaded all of the state and county information for US into the appropriate objects and created a point centered on the APSU campus. However, we can create a better map by filtering the data to focus on Tennessee and Montgomery County.

tn <- county %>% 
  filter(region=="tennessee")

montco <- county %>% 
  filter(region=="tennessee") %>% 
  filter(subregion=="montgomery")

With just these simple lines of code and included datasets we can use ggplot2 to quickly create a simple locator map for campus.

ggplot() + geom_polygon(data = state, aes(x=long, y = lat, group = group),
                        fill = "white", color="black") + 
           geom_polygon(data = tn, aes(x=long, y = lat, group = group),
                        fill = "gray", color="black") + 
           geom_polygon(data = montco, aes(x=long, y = lat, group = group),
                        fill = "red", color="black") + 
           geom_point(data = apsu_point, aes(x=x,y=y), color="black") +
  coord_fixed(xlim = c(-91, -81),  ylim = c(34, 37), ratio = 1.2) + 
  xlab("Longitude") + ylab("Latitude") + ggtitle("APSU, Montgomery Co., TN")

Keep in mind, that when using ggplot2 the order of the geoms is important for proper display of the image. Using the various options within ggplot2 you can customize the map based on your preferences.

2.1 Maps with raster images

It’s likely that the map you need to create is slightly more complex than the one above. Maybe you are interested in using some sort of imagery of basemap to add complexity or provide additional information for a locations. In the data folder of this exercise is a *.csv file called campus_points. This file contains the location of several points on the main campus. In order to view these files on aerial imagery of campus we are going to use OpenStreetMap to obtain a basemap for our closer look at campus.

Note: As of 2018, Google began requiring an API key for use with get_map that allowed for Google imagery to be used as basemaps. Although you can obtain a certain number of maps free per month, I have chosen to avoid that set-up for this example.

To use OpenStreetMap you need to first obtain upper-left and lower-right coordinates for the location you are interested in displaying. This can be done through programs like Google Earth or through min() and max() commands using the information in your dataset.

CampusMap <- openmap(c(36.5360,-87.3570),c(36.5300,-87.3495), type='bing')

APSU <- openproj(CampusMap, projection = "+proj=longlat +ellps=WGS84 +units=m +no_defs")

In the script above the bounding box was identified, the type of map listed (view ?openmap for additional options), and we set the projection for the map. Refer to the lecture for information rearding the projection and ellipsoid.

To add this imagery to the map we need to use the autoplot function in ggplot2. This function extends the plotting functions to include some spatial data classes. Otherwise the script for creating the map is similar to the simple map above.

campus_points <- read.csv(file = "./Data/campus_points.csv")
autoplot(APSU) +
  geom_point(data=campus_points, aes(x = X, y = Y, color=Name), size = 4, alpha = 0.8) +
  geom_text(data=campus_points,aes(x=X,y=Y,label=Name), color="black", vjust=-0.60, size=4.01, fontface="bold") +
  geom_text(data=campus_points,aes(x=X,y=Y,label=Name), color="white", vjust=-0.75, fontface="bold") +
  labs(x="Longtiude", y="Latitude") + theme(legend.position = "none")

With the addition of the aerial photo there is additional spatial information that can be obtained from viewing the map. Similar to the simple map above, you can customize the map through setting various themes and aesthetics in ggplot2.

2.2 Using KMLs or Shapefiles

Occassionally you might receive information from a collaborator that is not a *.csv file. Depending on the individual, that information might be provided as a *.kml file from Google Earth or in an ESRI shapefile. Using a simple import command in the rgdal package, these files can also be used with ggplot2.

To read either shapefiles or KML files you can use the following syntax:

readOGR(dsn="path to the shapefile",layer="name of the shapefile")

or

readOGR("path to kml/name of the file.kml")

With this in mind, you can easily use Google Earth to create your own datasets and import them into R for spatial analysis and mapping. This can be useful when collaborating with individuals or with citizen science projects.

In the data folder of this exercise you will find Campus_Points.kml file and a Main_Campus.kml file. The campus points file includes the same data as the previous *.csv file but as a kml, and the main campus file contains a polygon outline of the main portion of campus. You can use the script above to create datasets out of those files.

Working in ggplot2 with SpatialPointsDataFrame or SpatialPolygonDataFrame files can be difficult due to the structure of the data types. So to avoid adding additional packages or detailed explanation for this exercise, we will convert the points and polygon to a data.frame and spatialpolygon respectively. In later exercises we will work directly with these more difficult spatial classes.

To convert the point kml to a data.frame:

campusOGR <- readOGR("./Data/Campus_Points.kml")

OGR data source with driver: KML 
Source: "C:\Users\gentryc\Google Drive\APSU\Courses\BIOL 5700 Advanced Data Analytics\Exercise_9\Data\Campus_Points.kml", layer: "Campus"
with 23 features
It has 2 fields

campus_points <- cbind(campusOGR@data,campusOGR@coords)
campus_points[2] <- NULL
campus_points[4] <- NULL
colnames(campus_points) <- c("Name","X","Y")

This allowed us to import the kml, create a data.frame out of the coordinate and names, removed unnecessary columns, and renamed the columns.

To convert the polygon kml to a spatial polygon:

outlineOGR <- readOGR("./Data/Main_Campus.kml")

OGR data source with driver: KML 
Source: "C:\Users\gentryc\Google Drive\APSU\Courses\BIOL 5700 Advanced Data Analytics\Exercise_9\Data\Main_Campus.kml", layer: "Main_Campus.kml"
with 1 features
It has 2 fields

outline_points <- outlineOGR@polygons[[1]]@Polygons[[1]]@coords
colnames(outline_points) <- c("X","Y")
outline <- Polygon(outline_points)
sp_outline <- Polygons(list(outline),1)
outline_poly <- SpatialPolygons(list(sp_outline))
proj4string(outline_poly) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84")

This script created a spatial polygon and projected it to match the point data.

Now we can use these files to map the data from the kml files similar to the previous map of the csv data.

autoplot(APSU) +
  geom_polygon(data=outline_poly, aes(x=outline_poly@polygons[[1]]@Polygons[[1]]@coords[,1],
                                      y=outline_poly@polygons[[1]]@Polygons[[1]]@coords[,2]),
               alpha = .5, color="white") +
  geom_point(data=campus_points, aes(x = X, y = Y, color=Name), size = 4, alpha = 0.8) +
  geom_text(data=campus_points,aes(x=X,y=Y,label=Name), color="black", vjust=-0.60, size=4.01, fontface="bold") +
  geom_text(data=campus_points,aes(x=X,y=Y,label=Name), color="white", vjust=-0.75, fontface="bold") +
  labs(x="Longtiude", y="Latitude") + theme(legend.position = "none")

The script above is a little more confusing becase of the spatialpolygon class data. In order to map the aesthetics to the polygon we need to connect to the containers that hold the x and y values. This can be a confusing string, but we will discuss it further in upcoming exercises.

2.3 Student Submitted Mapping Question

How to create smooth polygons for species distribution maps from a set of points?

Following the “polygon kml to a spatial polygon” script above, install and load the smoothr package if you haven’t already. The smooth command in this package will take a polygon with geometric edges and “smooth” the data using one of the smoothing methods available in the package: chaikin, ksmooth, spline, or densify. Then you can adjust the smoothness value to determine how much smoothing is appropriate.

#Plot of the original polygon
plot(outline_poly)

#Smoothed version of the polygon
smooth_outline <- smooth(outline_poly, method="ksmooth", smoothness = 2)
plot(smooth_outline, col = alpha("grey", 0.2), border = "black", lwd = 1.5)

3 Interactive Maps

While there are a number of packages that can help create interactive maps such as plotly and mapview, for this example we are going to use leaflet. Leaflet allows you to create interactive web maps with the javaScript ‘Leaflet’ library. Because this package already has handlers for sf class data, we can use the original data we imported from the kml files to complete this portion.

In the most simple version, we set the point dataset we want to use for the map, tell leaflet the sort of markers we want to use to identify the points, and then let leaflet obtain the OSM tiles for the background image.

leaflet(campusOGR) %>% 
  addTiles() %>% 
  addMarkers(popup = campusOGR@data$Name)

Clicking on the markers will identify the place based on the “name” variable in our data. Because we have additional data from the polygon kml we can add to the map and customize it further.

leaflet(campusOGR) %>% 
  addTiles() %>%
  addPolygons(outlineOGR, 
              lng = outlineOGR@polygons[[1]]@Polygons[[1]]@coords[,1], 
              lat = outlineOGR@polygons[[1]]@Polygons[[1]]@coords[,2],
              fillColor = "grey",
              color = "black") %>%
  addCircleMarkers(popup = campusOGR@data$Name,
                   weight = 2,
                   color = "grey",
                   fillColor = "red",
                   fillOpacity = 0.7)

Here we customized the map with APSU colors and added the polygon to the background to identify the “main” campus area. We can take this a step further and provide additional information about the campus buildings and polygon by editing the attributes of each dataset.

outlineOGR@data$Description <- c("Austin Peay State University")

campusOGR@data$Description <- c("Academic", "Adminstrative", "Academic", "Library", "Academic", "Food", "Administrative", "Academic", "Residential", "Administrative", "Administrative", "Bookstore", "Academic", "Academic", "Residential", "Academic", "Academic", "Academic", "Academic", "Residential", "Residential", "Administrative", "Academic")

In the script above we provided a description of the polygon and descriptions of each building that can be used to identify the purpose of the building when the cursor hovers over the point.

leaflet(campusOGR) %>% 
  addTiles() %>%
  addPolygons(outlineOGR, 
              lng = outlineOGR@polygons[[1]]@Polygons[[1]]@coords[,1], 
              lat = outlineOGR@polygons[[1]]@Polygons[[1]]@coords[,2],
              label = outlineOGR@data$Description,
              fillColor = "grey",
              color = "black") %>%
  addCircleMarkers(popup = campusOGR@data$Name,
                   label = campusOGR@data$Description,
                   weight = 2,
                   color = "grey",
                   fillColor = "red",
                   fillOpacity = 0.7)

Notice now that when you move the cursor around there is additional information provided about each dataset. This can be used with numberous variable or to create an attribute table for each point when clicked.

3.1 Terrain and Elevation

Leaflet has a number of possible backgrounds that can be included as a basemap to your project. For example, if you were interested in elevation you can connect to a web-mapping service to obtain a zoomable terrain layer

leaflet() %>% 
  setView(lng = -112.956472, lat = 37.251343, zoom = 13) %>%
  addWMSTiles("https://basemap.nationalmap.gov/arcgis/services/USGSTopo/MapServer/WmsServer", layers = "0") %>%
  addMiniMap(zoomLevelOffset = -4) %>%
  addScaleBar()

3.2 Selectable Layers

Additionally, you can incorporate several basemaps and add an interactive selection tool to chose what background to display, set tranparencies for multiple backgrounds to be displayed at once, have live updates streamed to your map, or toggle the data off and on.

leaflet(campusOGR) %>% 
  addTiles(group = "OSM")%>%
  addProviderTiles(providers$CartoDB.Positron, group = "CartoDB") %>%
  #addProviderTiles(providers$Esri.NatGeoWorldMap) %>%
  addProviderTiles(providers$Esri.WorldImagery, group = "ESRI") %>%
  #addProviderTiles(providers$Stamen.Toner) %>%
    setView(lng = -87.353069, lat = 36.533654, zoom = 17) %>%
  addPolygons(outlineOGR, 
              lng = outlineOGR@polygons[[1]]@Polygons[[1]]@coords[,1], 
              lat = outlineOGR@polygons[[1]]@Polygons[[1]]@coords[,2],
              label = outlineOGR@data$Description,
              fillColor = "grey",
              color = "black") %>%
  addCircleMarkers(popup = campusOGR@data$Name,
                   label = campusOGR@data$Description,
                   group = "APSU",
                   weight = 2,
                   color = "grey",
                   fillColor = "red",
                   fillOpacity = 0.7) %>%
  addLayersControl(
    baseGroups = c("OSM", "ESRI", "CartoDB"),
    options = layersControlOptions(collapsed = FALSE),
    overlayGroups = "APSU")

This is only a very small fraction of the possible applications for mapping that can be done in R.

4 Your Turn!

Using the skills discussed in this exercise, create a site map for your thesis project or other dataset of interest. This map can be static or interactive, and can contain as many datasets or features as you feel are necessary to properly display your data. This map should be added to the project page you created in the last exercise to help build a “site description” section and/or “results” section if that data is currently available. Remember that this page should be consistently updated throughout the remainder of the semester culminating in your final presentation to the class in December. For an added challenge, you can try to fork this repository to give you access to the data files and a jumpstart on your script.